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COMPLEX EXPONENTIALS APPROXIMATION PROBLEM 
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Abstract. A kernel method is proposed to estimate the condensed density of the generalized 
eigenvalues of pencils of Hankel matrices whose elements have a joint noncentral Gaussian distribution 
with nonidentical covariance. These pencils arise when the complex exponentials approximation 
problem is considered in Gaussian noise. Several moments problems can be formulated in this 
framework and the estimation of the condensed density above is the main critical step for their 
solution. It is shown that the condensed density satisfies approximately a diffusion equation, which 
allows to estimate an optimal bandwidth. It is proved by simulation that good results can be obtained 
even when the signal-to-noise ratio is so small that other methods fail. 
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Introduction. Many difficult moments problems such as the trigonometric, the 
complex, the Hausdorff ones can be formulated as the complex exponentials approxi- 
mation problem (CEAP), which can be stated as follows, denoting random variables 
by bold characters: given a uniformly sampled signal made up of a linear combination 
of complex exponentials 
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where Cj , £j G (F, let us assume to know an even number n — 2p, p > p* of noisy 
samples 

d k = s k + e k , k = 0, 1,2, . ..,n- 1 

where e k is a complex Gaussian, zero mean, white noise, with finite known variance 
<7 2 . The CEAP problem consists in estimating p*,Cj,£j, j = l,...,p*. This is a well 
known ill-posed inverse problem often addressed in the literature, see e.g. [T3l [15l [TBI 

[H ig 1 m E E M ED M M [22] . 

We notice that, in the noiseless case and when p = p*, the parameters £j are the 
generalized eigenvalues of the pencil (U±,Uo) where U\ and Uq are Hankel matrices 
defined as 
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If we define Ui and Uo as Uq and U\ but starting from d^, k = 0, . . . , n — 1 it is 
evident that the generalized eigenvalues of the pencil (Ui,Uo) provide information 
about the location in the complex plane of the generalized eigenvalues £j , j = 1 . . . , p 
whose estimation is the most difficult part of CEAP. 
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To make precise this statement, let us consider apxp random matrix U and 
denote by ; j = 1, . . . ,p} its eigenvalues which form a set of exchangeable random 
variables. Their marginal density h(z), z e (Z*, also called condensed density [T5], is 
the expected value of the (random) normalized counting measure on the zeros of U 
i.e. 



If a pencil U = (Ui,Uo) of random matrices is considered, the condensed density 
of its generalized eigenvalues can be computed by the formula above. Its relative 
maxima provide information about the location in the complex plane of the generalized 
eigenvalues £j , j = 1 . . . , p. In [6l [5j [4] the use of this function for solving several 
moments problems was illustrated. However when the signal-to-noise ratio (SNR) 

c ■ I 2 

measured e.g. by SNR — miny is low the computation of the condensed density 
is very difficult even if we assume to have replicated observations as we do in the 
following. The main problem is that many relative maxima related to noise are 
present in the condensed density which are slightly less intense than those related to 
the true signals. As we are interested only on signal related relative maxima of the 
condensed density we look for an estimation method with noise filtering abilities. 

The aim of the paper is to show that it is possible to devise a kernel density 
estimation method which does have such properties. In [5] the same project was 
developed for the real exponentials approximation problem which is better known as 
the exponential analysis problem. The idea was to use a method proposed by Botev 
et al. [12], based on a class of kernels which satisfy a diffusion PDE, which allows 
to automatically estimate the optimal bandwidth. It was shown that the generalized 
eigenvalues of the pencil U can be approximated by the ratio of Gaussian variables 
whose density satisfies a specific diffusion PDE. In the limit for SNR I oo this PDE 
belongs to the class considered in [T3] and their method can be used. 

In the following a different approach is proposed. An explicit expression of 
the condensed density is derived and approximated by Laplace method. Then an 
anisotropic diffusion PDE is derived which admits the approximated condensed den- 
sity as a solution in the limit for a 0. Finally an optimal bandwidth is derived on 
the same lines of Botev et al. results. We can then propose an estimator which has 
better filtering abilities than the standard one based on Gaussian kernel. 

The paper is organized as follows. In Section 1 the condensed density of the 
generalized eigenvalues of U is derived. In Section 2 its Laplace approximation is 
computed. In Section 3 the anisotropic diffusion PDE is derived. In Section 4 the 
kernel estimator and the optimal bandwidth are computed. In Section 5 the proposed 
algorithm is illustrated. Finally in Section 6 two numerical examples are discussed. 
The Appendix contains the proof of most theorems and lemmas. 
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1. The condensed density of the generalized eigenvalues of U . From [51 
eq.(7)] the condensed density of the generalized eigenvalues of U is 



h n (z, er) = —E 
n 
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In the following Lemmas we give simpler forms and properties of the condensed den- 
sity. 

Lemma 1.1. 



{TT<7 Z ) n J a n/2-l Jan/2 



J* c = J£(C°'Ul), = n.((^,z), Qj = Q^Kz), Vj = v^\z\ 



and 



hk ~ 1 Jfc-1 



z , h = j 



X j (h,k) = x<ji ) k , qj = X j s, X,e(Z™/ 2x ", Q j = X j Xf eOJ n ^ Xn ^, 
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Lemma 1.2. If 
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X{h,k)=f h , s = X H c, v(Q = s H (I n -X H (XX H )- 1 X)s 
then u(Q = 0. 

Proof. = s H {I n -X H (XX H )- 1 X)s = c H X(I n -X H (XX H )- 1 X))X H c = 0. 



Lemma 1.3. When n — 2 the condensed density is given by 
h 2 (z, a) = e !+i*r 1 
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where p — denotes the SNR. 
Lemma 1.4. When n > 2 
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^ denotes the expectation with respect to the Gaussian density N (ji^Y,^ where 



and -y, /i are obtained by stacking the real and imaginary parts of 7 and [i. respec- 
tively. 

Lemma 1.5. 
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w/iere 0h(G^\z) are positive polynomials and Dj = det(Qj). Moreover 
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and P n /2{C^\ Z ) * s the only coefficient that does not depend on s. 
Corollary 1.6. When n > 2 
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where z and ( are the vectors obtained by stacking the real and imaginary parts of 
z and respectively, and 



1 ^ o ~,T\r< hM \tr~( h \ 8 I\ r ^\tr-~z\ 



-.18 



(T n 7r n/2 Z-^i 2 k 



k=0 



A- (C , *) 



Corollary 1.7. Mien n > 2, Vj = 1, . . . , n/2 

00) = /#>(*, 00) = 
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Moreover this is also the condensed density obtained when s — i.e. w/ien c = and 
it is circularly symmetric i.e. it depends only on \z\ 2 . 
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2. The Laplace approximation of the condensed density. The expression 
of the condensed density given in Corollary 11.61 can be approximated in the limit for 
(7 I as follows 

Theorem 2.1. In the limit for a 4 0, for z G Afj where Afj is a neighbor of ^ 
we have 

(2.1) hP{z,v) « h { i\z,a) = a- 2 G j (l U \z)e->^- U) ^ 

where 

G 3 (C U \i) = (2 7 r)f- 1 |^(C i , z)\-$Ki<l.,z) 



Hj is the Hessian of i/j and C is the unique minimum of Vj (£ , £j ) in a neighbor Dj 

~(j) 

°f {On h 7^ j} an d A)(C ! is given in Lemma \1.5\ 

3. The diffusion equation. In [6l Th.2] it was proved that h n (z,a) converges 
weakly to the positive measure ^ 2ft=i — 60 when cr J, 0. By the definition 
of hn\z,a) the same proof implies that hn\z,a) converges weakly to 5{z — 
Moreover from Corollary [O] we know that for a 4 oo, hn (z, cr) converges to a density 
independent of £. Therefore we can guess that when a moves from to oo, h„\z, a) 
diffuses from an atomic measure centered in £j to a measure circularly symmetric w.r. 
to zero. We then look for a diffusion equation which admits h„\z, a) as solution for 

0-4,0. 

Let us define z = x + iy and consider the anisotropic diffusion h t — L[h] where 



(3.1) L[h](x,y,t) =div 



h(x,y,t) 



and p(x,y) = hn\x,y, oo) is the stationary probability density, and aj(x,y) > is 
the unknown diffusion coefficient. Substituting h(x,y,t) in the equation above with 
hn\x,y,a), t = a 2 where now a is considered as a variable (not a fixed known 
value) and dropping the indices j, n and the variables (x, y) we get 

where 

Ei = -a y G y p 2 - a x G x p 2 + a y Gp y p + a x Gp x p + lap (G y p y + G x p x ) - 
a {G yy + G xx )p 2 + aG (p yy + p xx )p- 2aG {p 2 + p 2 x ) 



E% = a y Gp 2 v y + a x Gp 2 v x + 2aG y p 2 v y + 2aG x p 2 v x — 2aGp y pv y — 

2aGp x pv x + aGp 2 [y yy + v xx ) - Gp 3 
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In the limit for t 1 the dominant term in the left side of the equation is — 5 — 73-, 

p t 

therefore the equation is approximately satisfied when E 3 — or, equivalently when 
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> 0. 



In the following Lemma we prove that the Csiszar distance between h{x, y, t) and 
p(x, y) defined as 
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K x ^y^) 
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dxdy, $eC 2 :fi + ^ M+ , *"(■) > 0, *'(!) = 



is a monotonic decreasing function of t, therefore h(x, y, t) tends monotically to p(x, y) 
when t ^ 00 . 
Lemma 3.1. 



dt 



a{x,y)t>' 



p(x,y) 



p(x,y) 



dxdy < 0. 



4. The kernel estimator. Given a sample of size R of the data 

dW = [d«,...,dW], r=l,...,R 

where JE[dj c r ' ) ] = st, we consider the kernel estimator of hn (z,t) with bandwidth t 
given by 



1 - 



R 

r— 1 

where <&j(x,y,(,j(r)\t) is a solution of the diffusion equation (|3.1[) at time t with 
initial condition S(z — Cji r )) an d Cj( r )> J = 1, . . . ,n/2 are the generalized eigenvalues 
obtained from the data d^. Hence hn {x, y, t) is a solution of the diffusion equation 
(13. ip at time t with initial condition 

(4-1) E j (x,y) = -J2^-C j (r)) 

r—l 

Therefore Ej(x,y) is the empirical distribution of the generalized eigenvalue Cj and 

R 

^ P =i No- 
where /in ' (2) is the unknown true density corresponding to the known fixed value of 
a. 

In order to find the optimal bandwidth we need the form of the kernel <&j(x, y, d^; t) 
for t \. 0. By construction, dropping the index r, we have 

y, C; t) = y , t) = (x, y, C)e" * v * , U 0. 
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Associated to the anisotropic diffusion 13.11 there exists a Markov process X t whose 
transition probabilities are given, when t i 0, by &j(z, £; t). Moreover when the initial 
density of X is 5(z - (), the density of X t is ([H eq.(5.1)Ch.X.5]) 



Therefore the kernel $j(z, C; t)i and t > 0, satisfies the forward equation 

Moreover the conditional expectation of 8(u — X t ) on the hypothesis that Xo = z is 
([H eq.(4.5)Ch.X.4]) 
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and therefore the kernel $j(z, C; i) ; Vz and i > 0, satisfies the backward equation 
(4.3) 



^(^C;t)-L*[$i(«,C;«)] = o 
$ j (z,C;0) = 5(z-0 



where 
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and 



L[*j(«,C;t)] =div. 
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div,- 



V C Oi(C) 



is the adjoint operator of L and Ayv z and V z denote respectively the divergence and 
the gradient operators w.r. to the variable z. 

The mean integrated squared error (MISE) criterion to determine an optimal band- 
width t is given, dropping the index n, by 

MISE hU) (t) = E hU) [ [hW (z, t) - h& (z, t)] 2 dz 



'a 



MISE hU) (t)= [E h(j ){h^(z,t)}-h u \z)} 2 dz+ / Var h [h^(z,t)]dz. 



or 



(T 

Following Botev et al. [12] we have 
Theorem 4.1. 



[%){hW(z,t)} - h^(z)] 2 dz « i 2 ||L[/i«]|| 
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The MISE is then 



J Var h [h^(z,t)]dz « ^E h u [Gj\. 



MISE hU) (i) « i 2 ||L[^")]f + ^-E h u) [G 3 ] 



which has a unique real positive minimum in 
(4.4) tj = 



E h u) [Gj 



y AR\\L[h^)\\ 2 ' 
The optimal kernel estimator of hffl (z, t) is then given by 

^ ) (*.*i) = ^E*j( j8 'Ci( f )5*i) 

r— 1 

and the optimal kernel estimator of h n {z,t) is given by 

2 ™ /2 

(4.5) h„(z,ti,...,t n/2 ) - W^z,^) 

3=1 

5. The algorithm. In the following we assume that the generalized eigenvalues 
Q , j = 1 . . . n/2 of the pencils [u[ r \ Uq ] - where the Hankel matrices Uq , C/ 1 r ' ) 
are based on dY - have been computed for each r = 1, . . . , R and clustered in such a 

way that for each j = 1 . . . n/2 the j — th cluster is the set {(j r \ r = 1, . . . , i?} whose 
elements are independent realization of the r.v. £ •. The k-means method |17j can be 
used to solve the clustering problem. 

In order to use the optimal kernel estimator given in eqs. 14.41 and 14.51 for each 
j = 1 ... n/2 and r = 1, . . . , R we need to solve the initial value problem obtained 
by eq. 13.11 in the limit for a \. with initial condition given in eq. 14.11 By using 
the transformation h(x,y,t) — ~7~~|j the equation h t — L[h] can be rewritten as 

h f — -div 



aV(ji\ ■ Therefore the initial value problems are 



(5.1) 



where A is the Laplacian operator, Ej(x,y) is the empirical distribution of the gen- 
eralized eigenvalue £^ and 

a(x,y,j,r) = - P^v) v ^v) y jr { x ,y) =^{I n -Xf r (X jr Xf r )- 1 X j )s 

(t£(*.»)) + (t£(*.v)) 

where 

11 r=l { Z k , h = j 
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The density p( x, y) = hn (x, y, oo) is circularly symmetric (Cor. 11.71) . Moreover a 
closed form model of its modulus for each n can be found in [7] . 

The initial value problems were solved by a collocation method described in [2 3) 
in a m x x m y non-uniform grid. For each (J, r) the solution was approximated by 
the tensor product of Chebyshev polynomials in each spatial variable. Fast Fourier 
transform was used to compute the spatial derivatives. The resulting non-linear ODE 
system 

w'(t) = F(t,w(t)), w(0) = e, w(t),ee M m * m y 

was then solved by MATLAB's built-in function ode^.m, where F is the discretized 
right-hand side of eq. 15.11 and e is the discretization of Ej(x,y). The method is fast 
and stable provided that some spurious oscillations of the spatial derivatives close to 
the border of the integration region are filtered out. This task is accomplished by 
multiplying the derivatives by the function 

F(x, y) = [arctan{(x+7) /</>}— arctan{(i— j)/(j>}] [arctan{(i/+7)/0}— arctan{(y— 7)/^}] 
where 

X 7t{x ^mm)/ (,X max ^min) ^r/^, jj 7f(y Vmin) I (Vmax Vmin) 7f/2 

and 7 and <fi are suitable positive parameters. 

To compute the optimal bandwidth, E h u)[Gj\ is estimated by the sample mean of 
Gj(C, z) i.e. if the computed generalized eigenvalues are denoted by (j^ then 

1 R f 

®hi» [G 3 ] G 3 (Cl r \x, y)hg> (cj r) , x, y) dxdy a 

R m« m v ^■,.(c<'->,x h ,y fc ) 

i E E E &r(x h) y k , t 3 )f e *3 S x (h)6 y (k) 

r=l h=l k=l 

where 8 x (h) = x h - x h -i, S x (l) = S x (2), and S y (k) = y k - y k -i, S y (l) = S y {2). 
Moreover 



\\L[hP]\\ 



Q * 2 , R ra y 



dt 



^EEE 

r=l h=l fc=l 



d$ jr (x h ,yk,(j r \tj) 
dt 



8 x (h)5 y (k) 



To compute E h u)[Gj] and IILf/i^] || 2 we first need an estimate of tj which can be 
provided e.g. by the variance tj of the generalized eigenvalues in each cluster. 

6. Numerical results. In order to appreciate the advantages of the proposed 
kernel estimator, two numerical experiments were performed. R = 10 independent 

(r) (r) 

realizations d, ' = s k + el , k = 1, . . . , n, r = 1, . . . , R of the r.v. d*. were generated 
from the complex exponentials model with p* — 5 components given by 

p _ [ -0.1-i27r0.3 -0.05-i27r0. 28 -0.0001+i27r0.2 -0.0001+i2-7r0.21 -0.3-i2-7r0.35] 

c = [6, 3, 1, 1, 20] , n = 74, p = 37, a = 1, a = 3. 
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We notice that the frequencies of the 4 and 5 components are closer than the 
Nyquist frequency if n < 1/(0.21 — 0.20) = 100. Therefore a super-resolution problem 
has to be solved. To speed up the computations we limit the analysis of the condensed 
density to two regions containing respectively the first and second components and 
the third and fourth ones. The fifth component is isolated with a large amplitude 
therefore it is easy to identify even if its decay is fast. The considered regions are the 
rectangles defined by fl = tt x x fl y where fl x = [—0.8, 0.4], Vl y — [—1.4, —0.4] for the 
first region and fl x = [—0.1,0.6], Q, y — [0.5,1.3] for the second one. A mesh of size 
m x — 64, m y — 64 was considered. The values 7 = 1.6, <fr — 0-02 of the derivatives 
filter were selected by trials and errors. 

In order to apply the proposed method a pilot density estimate was first computed 
by the closed form approximation method given in [3] 



where R^ (z) is the R— factor of the QR factorization of the matrix U{ — zUq and 
/3 = 5na 2 [3j Prop. 6]. Then the generalized eigenvalues of the pencils [U[ r \ Uq], r = 
1, ...,R were pooled and the k-means method of clustering was applied with the 
number of clusters equal to the number of relative maxima of the pilot estimate. 

For comparison a Gaussian kernel estimate of density was also computed by Al- 
gorithm 1 in [T21 App.E]. 

In figure 16.11 the results obtained for the first region when a = 1 are plotted: the 
empirical density (top left), the pilot density computed by the closed form approxi- 
mation method (top right), the Gaussian kernel estimate (bottom left) and the result 
obtained by the proposed method (bottom right). In figure IB~2l the same results are 
plotted for the second region. In figures 16.31 and 16.41 the results obtained when a = 3 
are reported. The positions of the true complex exponentials are marked by a cross. 
It can be noticed that the proposed method is able to identify the two true complex 
exponentials even in the worst case, filtering out most of the spurious peaks of the 
empirical condensed density. Even if the location of the peaks is not perfect, it is the 
only method which is able to provide a reasonable solution to the super-resolution 
problem in the second region also for the smallest SNR considered. 
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Appendix. 

Proof of Lemma 1.1. Let us find /x. € @ n/2 , Q 3 € (F"/ 2xn / 2 and v } £ W such 

that 
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Choosing fj, ,Qj an Vj as defined in the text of the Lemma we get the result. 



Proof of Lemma 1.3. From Lemma [TTT1 we have 



where X\ 
therefore 



Qi 



l 



(ttcj 2 ) 2 
l + |z 2 | 



-i[(7-Ml) H Ql(7-Pl)+^l] 



Sq + ZSi 



l+\z 



d~; 



so 



si 



\s +zs!\- 

l+\z 2 \ 
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(w 2 ) 2 



/ | 7 | 2 e ^ d3?7d3 7 = 

'r 2 



(ttct 2 ) 2 l + |z| 2 V'^ 1 ' ' 1 + M ; 

If £ = — then 

i* = |so| 2 feC, and | m | 2 = | So | 2 |1 + " €|2 



l + |z| 2 ' ' Ul (l + \z\ 2 ) 2 

and 

ftWfc ff ) - e-^-^S ( N!il±M!. + I 

2t ' j " U^ 2 (i + kl 2 ) 3 ^(i + NI 2 ) 2 

But so = Ci and si = ci£i hence £i = £ and we get the thesis. 

Proof of Lemma 1.4. By considering the vector fi . obtained by stacking the 

real and imaginary parts of and \x , and the real isomorph Qj of the matrix Qj, 
and remembering that the Jacobian with respect to the real and imaginary part of a 
complex variable is Jr — |>/c| 2 , we get from Lemma 1 1.1 1 for n > 2 

h U) M = _}_f f |j*|2 e -*[(i-HP"Qi(i-Ap+^] dW o-)^ c o-) = 

(ttcH)" J Rn - 2 J Rn - - 

W /" n'bip n i^-&i 8 ni^-*i 8e "* [tt "^** a "^ 1< ^ Wd9 c w - 



(trt 2 ) 
By defining 



l,n/2 



and noticing that 

| 7l | 2 =7 T A,7 

we have 

JR™- 2 {nv ) n/ \Qj{C , z )\ 2 

and the thesis follows. 

_ i 

Proof of Lemma 1.5. By considering the normalized vectors y. = 2 7, we 
have y ~ AT (rtij , i) where 
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where Dj — det (Qj) and Qj — adj(Qj). Moreover 



where 



Bi = Sf^sJ = ^-Q7 1/2 A i Q7 1/2 



^q) /2 a iQ -) /2 

'i 



From 1, Th.l], denoting by Qi the quadratic form yjBiy^., we have the recursion 



l,n/2 

1=1 



i/2-l n/2 



,/2 



1=0 j 1= 2 j i= 2 



Q 2 • • ■ Q„/ 2 



2ji ■ • ■ 2.. 



where for i = 0, g = m T Bim + tr(Bi) = E[Qi] and for i > 0, ji ^ j 2 ^ ■ ■ • ^ ji and 
9n-k = m T {BxB h . . . IS, ■ IS, I!. /,',. ...B k +,. .+B h B h . . . B Ji B 1 )m+tr(B 1 B n ...B u ). 
But then 



9ji. 



Tlj A lj + jnji tr ( A iQj A hQi ■ ■ ■ A kQi) 



T 2( 



531. 



= 217 [^d q1 j' 2 S a 1 + H A iQj A 3iQj ■ ■ ■ A nQ.i)^ 

where 

A = O l .\ ] Q J .\,Q, . . . Q,.\, (), ■ (),A, Q,.\:Q,.\,.Qj . . . A ji Q j +. . .+Q j A n Q J . . . A ji Q j A 1 
We have 

3 3 

E[QiQ 2 ] = #[Qi]£[Q 2 ] + 4mjB 1 B 2 m j + 2tr{B 1 B 2 ) = 



JytgJ Qj A l Qjgjfj Qi A 2 Qjg 



2 

^3 (^Q j A 1 Q j q j tr{A 2 Q j ) +^Q j A 2 Q j q j tr{A 1 Q j )+^Q j A i Q j A i Q j g :j 



AD* 



tr{A x Q j )tr{A 2 Q j ) + tr(AxQjA 2 Qj) 



and in general 



«/ 2 Jlh 
fj(C«\z,a) = P n/2 (a\z) = Y J h{&\z)—— K 



2 h D r ' 

where (3h(C^\z) are positive polynomials. Moreover P n / 2 (C,^\ z) is the only coeffi- 
cient that does not depend on q and therefore it does not depend on s. 
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Proof of Corollary 1.6. By Lemma H~4l 

- (TTa 2 ) n / 2 \Q.j{Q j) ,z)\2 



2 fe J R n-2 ~ ' n "-fc+i - - 



fc=0 

SI8 



Noting that |£ r - Chi = ICr - Cftl and iCr - A = \Q r - z| 8 , we get the thesis. 

Proof of Corollary 1.7. Noticing that in the definition of gj(( ,z,a) when 
a 1 oo all terms vanish but the last one, we get the first part of the thesis. By Lemma 

11.51 8„f)(C ,z) is the only coefficient which does not depend on s. Therefore when 
s = equation 11.11 reduces to hn (z, oo). Finally by symmetry the condensed density 
does not depend on j, therefore all hn\z,co) must be equal. Moreover in [7] it was 
proved that, when s = 0, it depends only on \z\ 2 . 

Proof of Theorem 2.1. We recall that if 

1=1 K{y)e- Xu{ y ] dy, V open set C M d , A € M + 

Jy£V 

and v(y) has a unique minimum in T> and this minimum occurs at a stationary point 
y of v(y), then the Laplace's approximation to I is given by 

1 = {2ir)i\-$\H{y)\-^K{y)e- X " { y) 
where H(y) is the Hessian of v. 

We know that Vj > and, by Lemma fl.2[ £ is the only vector such that Vj{^) = 0. 

~U) ~ . . * ~~ 

Therefore by continuity, has a unique minimum £ in a neighbor 2?j of 

{0i> ^ 7^ i}- Moreover from Corollary 11.61 we notice that the dominant term in the 

~(J) 

sum defining <?j(C , 2, <r) when a 1 is the first one, therefore in this case 

/?&*)- \ 1 o /?CJ) ^nr^VilCr-Ca^rir-^-ICr-^l 8 1 -(j) 

■Dj (C j«) 

Then by using Laplace's approximation with A = and (i = n — 2, we have, for 
z e A/} where A/j is a neighbor of ^ 

Di 



^■)(i,a) = 4r / #i(C (i) , 



(6.1) « C 7- 2 (2 7 r)f- 1 |^(e.^)|-^ j (C i ,5)e" >i(£ ^" ) 

where Hj is the Hessian of Vj. For simplicity we will denote this approximation by 
(z, a). Let us define 

G,(C U \S) = (2n^- 1 \H J (C,z)\-iK 3 (C,5) 



the same symbol hn (z,a). Let us define 



then 

(6.2) hW{z,a) = a- 2 G J (^\~z)e-^^ ) ' i \ 
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Proof of Lemma 3.1. 

&D(h,P) 
dt 



R 2 at 



h(x,y,t) 

p(x,y) 



*' (^4-)h t (x,y)dxdy 



a(x, y)V 



h(x,y,t) 

p(x,y) 



dxdy. 



Integrating by parts we get 
dV(h,p) 



dt 



( div *' ( 

Jr 2 l V 



*' (^4] a{x,y)W ( h{X ^ t] 



P(x,y) 



P(x,y) 



dxdy- 



p(x,y) 



— ,' h(x,y,tY . 
■a{x,y)V\ ^ )dxdy, 



where • denotes scalar product. By the divergence theorem the first integral is zero 
because ^'(l) = 0. Moreover 



p(x,y) 

Proof of Theorem 4.1. 



, ( h{x,y,t) 



„ ^ h(x,y,i) ^ v ( h(x, y, t) 



p{x,y) 



p(x,y) 



By definition of adjoint operator, taking into account that lirn^oo $j(z, (; t) = we 
get for t 4 



dE hW {bW(z,t)} 



dt ja 



L*[<s> j (z,(;t)]hW(0d<;= [ ^(z,C;t)L[h^(0]dC 

Ja 



But lim t ^o ®j( z i Ci t) — S(z — () hence 

/ ^(z,C;t)L[h^(Q]dC ^ L[h^(z)}, no. 
Ja 

By considering the first order Taylor series approximation of E h a){h^(z,t)} we get 
E W) {h« (z,t)} = E hW {h« (z,0)} + 1 ^ilMlMl | t=0 + (t% 



but 



hence 



and 



E hW {hW(z,0)}= [ $ j (z,C,0)h( j \()d{= f 5(z-QhV\Qdt = hW(z) 
Ja Ja 

E hU) {h«\z,t)} = h^(z)+tL[h^(z)}+0(t 2 ) 

[ [E h(j) {h^(z,t)}-h^(z)] 2 dz^ [ [tL[h^(z)]] 2 dz^t 2 \\L[h^]\\ 2 
Ja Ja 
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For approximating the integrated variance let us consider first the second moment 

E hU) {$ j (z ) C,t) 2 }= f <s> 3 (z,(;t) 2 h^(()dc= f r 2 G^(*,0VH<*-OfcW(C)dC 
Ja Ja 

But it was proved in [6j Th.2] that 



t-yO t 



hence, for t J, 0, 



m 



As C-i( r ) are independent Vr., it follows that 



Kor h [h^ (2, t)] =Var h 



^ i? l fl 



r=l 



±E m [<bj (x, y, C,- (r); t) 2 ] - | [$, (*, j,, ^ (r); t)]) 2 



^G,(z,z)^(z) - (^)(z) +^[^)(z)]) 2 « ±- Gj (z,z)h^(z) 
because for t \. the second term is negligible w.r. to to the first one and 
Jvar h ^{z,t)\dz» ■±jG j {z,z)hU)(z)dz = ±E hU )[G j ]. 



17 





Fig. 6.1. First region, a = 1. The empirical density (top left), the pilot density computed by 
the closed form approximation method (top right), the Gaussian kernel estimate (bottom left) and 
the result obtained by the proposed method (bottom right) 




Fig. 6.2. Second region, a = 1. The empirical density (top left), the pilot density computed by 
the closed form approximation method (top right), the Gaussian kernel estimate (bottom left) and 
the result obtained by the proposed method (bottom right) 
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Fig. 6.3. First region, a = 3. The empirical density (top left), the pilot density computed by 
the closed form approximation method (top right), the Gaussian kernel estimate (bottom left) and 
the result obtained by the proposed method (bottom right) 




Fig. 6.4. Second region, a = 3. The empirical density (top left), the pilot density computed by 
the closed form approximation method (top right), the Gaussian kernel estimate (bottom left) and 
the result obtained by the proposed method (bottom right) 



